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ORIGINAL FAGK IS 
OK HOOK QUALITY 

1. ic;rcdugLisa a a a aiaiarlaala 


The object of the present paper is to introduce an expli- 
cit integration scheme for the fculer equations of gas dynamics 
which, after extensive tests, 1 consider the best among all 
schemes available so far. 

According to P. Cordon, from whom 1 got the first sugges- 
tion of the scheme UJ, its basic idea could be found In an early 
paper by Courant, Isaacson and Pees [2J. In 19b7, after doing my 
homework on Cordon's scheme [3J, I dropped it, giving my prefer- 
ence to MacCormack's scheme 14] which was also welcomed by a 
large number of numerical gas dynamicists. 

Three reasons prompted me to adopt MacCormack's scheme 
instead of Coraon's: 

1) The former has second order accuracy, the latter had 
only first order accuracy, 

2 ) The former could be easily coded for any number of 
space-like variables, the latter was limited to one-dimensional 
problems, and, on a minor scale, 

3) The former seemed to need simpler codes. 

A couple of years ago, L. Zannetti rediscovered tne basic 
idea of the integration scheme independently. He was, however, 
in a more advantageous position than Goraon because he had a good 
experience of two-level schemes and he knew how to maintain 
second order accuracy despite using one-siaed approximations to 
space derivatives [ 5 j . The first public appearance of Zannetti’s 
scheme lb] snowed that the first objection above had been re- 
moved; the scheme, however, was still essentially one- 
dimensional. 

Since then, 1 have been working on the scheme, in order 
to test it in the widest possible range of flow parameters ana 


- 1 - 



Introduction and historicals 


under the most exacting circumstances. I did also work out a 
multi-dimensional version of the scheme and 1 consider, thus, 
that the second objection above can also be removed. As far as 
the third objection is concerned, the code turned out to be not 
as complicated as I thought and it can be provided now in the 
form of a 'black box' which maxes its application very simple. 

Therefore, 1 believe that the scheme whlch’l am going to 
describe deserves recognition. It is superior to any other from 
a physical standpoint and it provides results of great accuracy. 
The second statement is an obvious consequence of the first! 

To identify the scheme, 1 decided to call it the 
X-scheme, since 1 could not attribute to it a well-defined pater- 
nity. In fact, most of our recent progress in numerical analysis 
seems to be the product of a natural evolution, occurring in dif- 
ferent places at almost the same pace. A paper by F. Walkden, P. 
Caine and G. T. Law;. . which Just appeared in print [17], seems to 
overlap a good part o ** the present one although the Authors may 
have missed what I consider a very important improvement over the 
Gordon scheme. The scheme, as outlined by halkden, has indeed 
only first order accuracy, and seems to make heavy use of averag- 
ings whose damaging influence cannot be evaluates « priori but 
wnicn smell strongly of artificial viscosity. 


2- Recasting &X lule.ii ' a equations 


In what follows, we will denote by T a time-like coor- 
dinate, and by X, X two space-like coordinates. Such coordi- 
nates are in general not Cartesian, but obtained from Cartesian 
coordinates through a series of mappings and stretchings, as re- 
quired by the problem. 

As the main thermodynamical unknown we will assume the 
logarithm of pressure, P. In time-dependent problems, we will 
assume, as the main kinematic unknowns, the velocity components 
along the computational grid lines in the physical plane. In 
steady, supersonic flow problems, it is more convenient to assume 
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ORIGINAL PAGK IS 
OF POOH QUALITY 


the slope of the velocity vector or two angles defining the velo- 
city vector. 

beginning with one-d iaens ional problems, that is, prob- 
lems whicn contain only one space-like variable, the typical un- 
steady flow is defined by the equations: 


P ♦ a P ♦a u ♦ c s 0 
T 11 X 12 X 1 


u ♦a P ♦a u ♦ c s 0 
T 21 X 22 A 2 


and tne typical steady, supersonic fLow is defined by similar 
equations, where u plays the role of the velocity slope. The 
coe f f ic ients , a depend on the computational grid. Tne terms, c 

may also depend on other conditions or assumptions; for example^ 
in unsteady flows c^ is not zero if cnanges in a duct cross- 
sectional area are considered, c. is not zero if the frame of 
reference is accelerating. 


The simplest possible case of insteady flow is obtained 
if the frame or reference is fixed and Cartesian, if there are no 
area effects and no stretching of coordinates is used. In that 
case, for a perfect gas with a constant ratio of specific he. ts, 
y , M) become 

P ♦ UP ♦ YU * 0 
T X X 

( 2 ) 

2 
a 

u ♦ P ♦ uu s 0 
T Y X X 

Similarly, the simplest possible case of steady, supersonic flow 
is obtained if the frame of reference is fixed and Cartesian and 
no stretching of coordinates is used. In this case, if (for 
better clarity) we write a instead of u for the slope of the 
velocity vector, and T and X are the horizontal and tne vert- 
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ioal Cartesian axes, respectively, (1) become 

P ♦ (o/x)P ♦ (y/c)o » 0 

2 2 2 ( ^ ) 
a ♦ (a /yu )(1 ♦ o /c)P^ ♦ (e/c)e^. » 0 

where u is the velocity component in the Indirection, and 

k « 1-a 2 /u 2 (4) 


The slope of the characteristics of system (1) in the 
(X,X) plane, X s dX/dT, is obtained by solving the equation: 





(5) 


There are, in general, t*.o roots to this equation. he will 
denote by X^ and X^ the smaller and the larger root, respective- 
ly. It is well-known that the system (1) car. be replaced by the 
equivalent system: 


( a -X 

)(P „«.X P ) 

- a ( u ♦ X u ) 

♦ (a„ -X )c. 

-a _C_ 

M 

o 

22 1 

T 1 X 

12 T 1 X 

22 1 1 

12 2 


( a„ -x„ 

)(** *X P ) 

- a . ( u u ) 

♦ ( » -X ) c . 

- a e 

S 0 

22 2 

T 2 X 

12 T 2 X 

22 2 1 

12 2 



Obviously, if tne second equation (6) is subtracted from the first 
and the second equation, multiplied by a -X , is subtracted 
from the first, multiplied by a -X • we obtain Cl) again. If, 
however, we write P and u^ for P^ and u^ in the first of (o) 
and, similarly, P^ ana u in the second, we obtain 



A ) P -A X P 
1 1 XI 2 2 X2 


D (X P -X P 
21 2 X2 1 XI 


+D (Xu -Xu )+c 
12 2 X2 1 XI 1 


♦ A X u -A Xu ♦ c 
1 2 X2 2 1 XI 2 


= 0 
X 0 


(7) 


where 
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lo ) 


trci the viewpoint of partial differential equations, P a p 

X 1 X2 

and u^, * u^ ^ ; it is easy to verify that (7) are, once more, 

identical with (1). Not so from tne viewpoint of finite differ- 
ence equations, if we decide to approximate P ana P in dif- 
ferent ways, and similarly do with u^ and u^* 


Again, in the simple unsteady flow defined by (2), 


a n " u • 


*12 * ^ ’ 


* u-a , 


*2 1 * 8 /T 1 


X^ * u*a 


a a u 

22 


(9) 


and (7) become 


P ♦ ( 1 P 

T 2 1 XI 


X„P V ,) V T~ (X u - — Xu v , ) = 0 
c ke. «.a 2 x2 1 XI 


( 10 ) 


U ♦ (X p -X P ) ♦ (X u ♦Xu ) r 0 
T 2 Y 2 X2 1 XI 2 2 X2 1 XI 


whereas, in the simple supersonic flow defined by (3), 


\2 

1 8 

a s <3 /< , a S Y/< , a a , a a O/K 

11 ’12 ’21 4 22 


YU * 


1,. a-A, 


. a-JL, 


! ..2 


X. a -(o-“) , X, a -‘■(0 ♦ “) , 0 a \ ! M~ - 1 


1 K ' 2 

U 


2 K 


(11) 


and ( 7 ) become 


P ♦ ^(X P ♦ X p )- 4_a (I o -X a ) a 0 
I 2 1 XI 2 X2 2 2 X2 1 XI 

2 YU 


( 12 ) 


o ♦ 1 * ( X P -X P ) ♦ - L (l o ♦ X o ) a 0 

T , 2, 2 (2 1 XI 2 2 X2 1X1 

2a 6 
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i. flla.srftU&aLlan 


from t n * two examples abo'e, it is clear that 1^ and 1 
are the slopes of the two chara x eristics which carry the it for- 
mation necessary to define P and u (or P and a) at every 
point. For auch a point, A (Fig. 1) the domain of dependence is 
defined by point B and point C, the two characteristics being 
defined by AC and AB, respectively, to a first degree of accu- 



racy. When discretizing the partial differential equations for 
computational purposes, point A must be made depend on poi'.ts 
distributed on a segment which brackets BC [7], for example on 
points D,E and F of Fig. 1. Such a condition is necessary for 
stability but it must be interpreted with a grain of salt. Sup- 
pose, indeed, that one uses a scheme in which a point such as A 
is always made depend on D,E and F, indiscriminately (this is 
what happens in most of the schemes currently used, Including the 
MacCormack scheme). Suppose, now, that the physical domain of 
dependence of A is the segment BC of Fig. 2. The information 
carried to A from F is not only unnecessary; it is undue and 
therefore the numerical scheme, whilst not violating tne Cr'L 
rule, violates a law of physics, the law of forbidaen signals. 
Physically, it would be much better to use only information from 
0 and E to define A, even if this implied lowering the nomi- 
nal degree of accuracy of the scheme. In other words, to say 
that a given scheme, using points 0, E and F, has a second 
order accuracy is meaningless since a wrong scheme has no accura- 
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cy whatsoever. 


or ' K '' 




In transonic problems treated by relaxation letnoas, the 
need for a switching in tne discretization scheme in passing froa 
subsonic to supersonic points was felt from the beginning Co], 
and it is evidently related to the law of forbidden signals, 
altnougn we find it Justified in the literature as a gimalcx to 
avoid expansion shocxs or to give the proper C rection to the in- 
crease in entropy [9] (a curious stateaent, indeed, for a tech- 
nique which is, by definition, isentropic ! ) . 



Fig. 2 

The sensitivity of results to the nuaerlcal doaain of 
dependence as related to the physical doaain of dependence ex- 
plains why computations which use integration schemes sucn as 
hfiCormack's show a progressive deterioration as the AC line of 
Fig. 1 tends to become parallel to the T-axis ( x ^ »0 ) , even if x ^ 
is still negative [10]. The information frca F actually does 

not reach A; in a coarse aesh, such information may be drasti- 

cally different froa the actual values (froa C) which affect 
A. On the other hand, since the CFL rule must be satisfied and 
F is the nearest point to C on its right, the weight of such 
information should be minimized. The X-scheme, relying on (7), 
provides us with such a possibility. 

We stipulate that every derivative f be approximated 

by using points which lie on the same side of £ as C, and 

that every derivative f be approximated by using points 

which lie on the same side of £ as B. by so doing, not only 

we relate each characteristic witn information which is found on 


7 



Discretization 


tn« same aide of A from which the characteristic proceeds, but 
we weigh such inform, ton with factors whloh contain 1 ^ and so 
that the contribution of points located too far outside the phy- 
sical dcaaln of dependence is minimized. 


A one-level scheoe which defines 


( f -f .)/AX 
F Ik 

(f ,-f)/AX 
£ D 


(^< 0 ) 


(^> 0 ) 


(U) 


is Gordon's scheme [1], accurate to the first order. To obtain a 
scheme with second-order accuracy, we consider two levels, in a 
manner very similar to MacCormack ' s . hie must, however, introduce 
more points, as 'n Fig. 3* At the predictor level, we define 


(f F‘ f E )/AX <V°> 

f Xi * 

(2f -3f. ♦ f _ ) / A X (1 >0) 
a D u 1 

At the corrector level, we define 


(14) 


<-2f ♦ 3f -f D )/AX ( A <0 ) 

A N P X 

f xi • 

(f -f)/AX U>0) 

AM X 

It is easy to see that, if any function f is updated as 


(15) 


? » f ♦ f AT (16) 

T 

at the predictor level, with the T-derivatlves defined as in (7) 
and the X-derivatlves defined as in (14), and as 

f ( T ♦ AT) * -f(f ♦ T ♦ f AT) (17) 

. 2 T 
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ORIGINAL Ra>. 
OK P(X)R UUAI.IJ 


at tne corrector level, witn the T-der ivatlves aefinec again as 
in (7) ana the X-aer ivat ives aefinea as In (15), tne value of f 
at T ♦ AT is oDtainea with secona order accuracy. The updating 

T 1 

M A N P 

O o o o o 

AT 


-o 

G 


o 

D 


o- 

E 

Fig. 3 


o 

F 


X 


rules (16) ana (17) are the same as in the MacCoraack scheme. 


1. Addition *1 remarks un A ' blac k Bax' mail 


In the simple examples above, the computational grid 
coincides with the physical grid, which is fixed. The charac- 
teristics are the physical characteristics of the flow. For ex- 
ample, in the unsteady case, the domains of dependence are relat- 
ed to the simple concept of subsonic flow (for u>0, u<a, 
A <0.X >0) and supersonic flow lu>a, 1^>0,1 >0); similarly, in 
the case u<0. The system (1), however, ar well as (7), is ouch 
more general. The computational grid, in the (X, T) plane, is 
fixed but its counterpart in the physical plane generally varies 
witn T. Tne discussion of the domain of dependence is carried 
on and the slope of the characteristics is determined in the 
(X,T) plane, with reference to the fixed grid. Therefore, it may 
happen that a flow which is physically supersonic has to be con- 
sidered as 'subsonic' in the computation (Fig. 4); conversely, a 
physically subsonic flow may have to be considered as 'superson- 
ic* in the computation (Fig. 5). 

We conclude our analysis of the i-scheme by showing (in 
Fig. 6) an example of a FORTRAN subroutine to compute the X- 
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Fig . U 




derivatives according to the conventions outlined above. There, 
LOOP*- 1 , 1 Ln the predictor and corrector stage, respectively, 
N cenotes tne grid point location, LAPs(L00P ♦ 1)/2, and 
ALAMO), ALAM ( 2 ) are -X 1 and X^, respectively. Special provi- 
sions are made at Ns1,2 and at NsNA, NC where NC is the last 
grid point and NAsNC-1 to account, in this case, for the pres- 
ence of two rigid walls. 


1- lR£ BQUS.Ll.Qa SLl ZA.&C&1 


The equation of energy is conspicuously absent froo our 
analysis, so far. Indeed, written in terms of entropy, it simply 
states the Invariance of S on a particle path: 
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The equation of energy 


UKiiiJNAI, 

<* **>« guAiSw 


X-l 

NTL-N-LOOP 

LX>LAP«1 

LX1-3-LX 

no-n-loop 

Nl-N*LOOP 
N2-N1+LOOP 
PN-P (N) 

UN -U (N) 

PNO-P(NO) 

UNO-U (NU ) 

PNl-P(Nl) 

UN1-U (Nl) 

PN2-P(N2) 

UN2-U (N2) 

IF (NTL.NE.-2.AND.NTL.NE.NAICO TO 9 
PN2-3.* IPN1-PN)>PN0 

« 


UN2* J . • (UN 1 -UN I ♦L'NJ 
9 IF (ALAN ( I.X I . LT. t )GO TO 1 

IFINTL.EC. 1.0R.NrL.Ew.-SC)G0 TC 3 
PX(LX)- IPU-PNO) *DDXD 
UX(LX) • (L'N-UNu) *DOXU 
IF(NTL.EQ.-1.0R.NTL.EQ.NC)G0 TO 7 
I F ( ALAM ( LX 1 1 . LT . 0 . ) GO 70 7 
LSAV-LX 
LX-I.X1 
LX1-LSAV 
1 I "LX ♦ LAP 

3 PX(LX) • (-2. *PM*3.*PN1-PN2) *DDXD 
UX(LX)- <-2.*UN»3.*UNl-UN2) *DOXO 
IF (I .tO. 2) RETURN 
7 PX(LXl)-PXUX) 

UX (LX1 ) -UX (LX) 

RETURN 


Fig. b 

DS/Dt*0 (16) 

where t is real time (In the two simple examples above, (1o) 
becomes 

S„. ♦ uS «0 (19) 

T X 

and 

S ♦ os si) (20) 

T X 


respectively) . 


At each integration step, the other equations of motion 
can and should be integrated independently of (16); the latter, 
in turn, should be integrated by approximating the X-derivative 
of 6 in such a way that information is never carried backwards 
on a particle path [11]. Obviously, the rules for the approxima- 
tion of 6^ are different from the rules for the approximation 



Tne equation of energy 


of u^, and any conceptual or coding mixup should be care- 

fully ’ avoided . The only interaction between entropy and the oth- 

occurs via the value of 
,* in steady flows (see 

(3)) since, in general, 

a 2 » y exp[(t-1)P/Y ♦ S/y] (21) 


er dependent variables of a flow probl 
a" in unsteady flows (see(2)) or a 


em 


and, for steady flows, 


q 


2 


JL-. 2 2 

, (a -a 
t - 1 o 


) 


(22 ) 


where a is the stagnation speed of sound, and q is the modulus 
of the velocity. 


fcxlnumm is. aul-l.lfl.ungQ a legal problems . Unsteady flows 


Let us consider now a problem with two space-like vari- 
ables, X and Y. The typical unsteady flow is describea by the 
equations : 


P T 

* a >, p x 

♦ 

*!2 U X * b 11 P Y 

♦ 

b 12 V Y 

♦ 

g u 
*11 X 

♦ 

X 

> 

<M 

to 

* C 1 

* 0 

U T 

* *21 P X 

♦ 

a 22 U X * “ 2 1 P I 

♦ 

h 22 V Y 

♦ 

Vy 

♦ 

£ V 

*22 X 

+ C 2 

* 0 (23) 

V T 

* B 31 P X 

4- 

"32 U X * b 2l P J 

♦ 

b „ „ v 
22 1 

♦ 

S i. U Y 

♦ 

d v 
3 X 

# °3 

« 0 


The 

simplest possiole 

case 

Of 

unsteady 

flow 

is 

obtained 


if the gas is perfect, the frame of reference is fixed and Carte- 
sian, and no stretching of coordinates is used. In that case, 
(23) become 


P T * UP X * VP Y * TU X * Y V° 

2 

« T ♦ uu x ♦ vUj ♦ (a /y)P x »0 


V T * UV X * VV Y * (a /y)P y*° 


( 24 ) 
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OKkUNAI, I’AliK lh 
OF P»K)H QUALITY 

be see that not only the terms Dut also the terms affected by 
the coefficients g and n depend on the choice of the compu- 
tational g-_d. for a well-chosen grid, we expect such terms to 
play a minor role in the computation. Host of the remaining 
terms can be arranged in two matrices, one formed with X- 
..erivatlves of P and u, the other with Y-derlvatlves of P and v . 
Two terms are left out, d,u^ and d^v^, which, as seen from (24), 
are physically relevant. 


be will now rearrange the equations (23) in two sets, one 
essentially related to X and the other essentially related to 
I: 


and 


P ♦ a P 

T 11 X 


a u ♦ * u ♦ c * 0 
12 X *11 Y 1 


“t * *J1 P * ♦ *2J U X * V* * V* * V° 


(25) 


P ♦ b P ♦ b v ♦ g v s 0 
T 11 Y 12 Y S 1 2 X 


v ♦ b P ♦ b v ♦ d v ♦ h P ♦ c *0 
T 21 Y 22 Y 3 X 3 X 3 


(2b) 


X Y 

The meaning of P^ and P is simply that P . in (23), is intended 

T i T X 

to be obtained as the sum of two separate contributions, P and 

P : cherefore, if the first of (25) and the first of (2b) are ad- 
T 

def together, one obtains the first of (23) again. The distribu- 
tion of the terms, and 0 A * between the first of (25) 

and the first of (2b) is arbitrary, but the arbitrar iness is ir- 
relevant, as we will see soon. 


Both (25) and (2b) show 
ticularly if we rewrite them in 


a strong resemblance to (1), 
the form 


par- 




P„ 

♦ a u 

♦ C * 

1 1 

X 

12 X 

1 


P 

♦ a u 

♦ c x . 

21 

X 

22 X 

2 


0 

0 


(27 ) 
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and 


with 


(26) , 


and 


where 


P J * a u P r * Sj’i ♦ C I*° 

Y 

v ♦ b P ♦ b v ♦ c » 0 

T 2 1 x 22 y 2 


(2b) 


„X 

c .*‘n“> * c i • 

c i**i2 v x 


C i* a 2 U « * V» * C 2 


W* * V* * e 3 


(29 ) 


Tnerefora, we can proceed as in Section 2 for (27) and 
separately, to obtain two sets of equations similar to (7): 


XX X XXX X 

P ♦ A x P -A x P ♦ D (X u -x u ) ♦ c *0 
T 1 1 XI 2 2 X2 12 2 X2 1 XI 1 

u + 0 (IP -IP ) + A X u -A Xu ♦ C *0 
T 21 2 X2 1 XI 12 X2 2 1 XI 2 


(30) 


y y Y 

P ♦ B X*p -B X p 
T 1 1 Y1 2 2 Y2 




y y y Y Y Y 

v ♦ D (X P -X p ) ♦ B X v 2-B X v 1 ♦ C «0 
T 21 2 Y2 1 Y1 1 2 y 2 1 Y 2 


(31 ) 


. X ,Y 

a -X p -X 

. -LI i _ 111 i 

As , D s 

i X X i 


X X 
X -1 
2 1 


x X -x Y 

2 1 


_X a i1 Y - b i.1 

D iJ * X .X ’ ij * 


X -X 
2 1 


y y 
x -x 
2 1 


(32) 
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X J 

and A^,A^ ar * solutions of equations similar to (5): 


a -X X 

a 


Y 

b -A b 

1 1 i 

21 


111 21 

a 

.X 

a -A 

«0 , 

b b -A* 

12 

22 i 


12 22 1 


>0 


( 33 ) 


In discretizing (30) and (31). we will use the same rules as 
stated in Section 4. 

he turn our attention now to the discretization of the X- 

X 

and Y-der^vatlves in (29). The most important terms, u in C 
and v^in C ^ , should be discretized recalling that they come from 
Lagranglan derivatives in momentum equations. Therefore, it 
seems proper to discretize them the same way the entropy deriva- 
tives are approximated in the energy equation. Specifically, the 
sign of d ; must be tested; then, at the predictor level, 


V s 

X 


(V -vJ/U 
F E 


( 2 V 3 W /AX 


(d 3 <0) 


( d _ >0 ) 


(3*) 


and, at the corrector level, 


( - 2v « * 3 vV'“ 


(v -v ) / A X 

A M 


( d ^<0 ) 


(d >0 ) 


(35) 


Similar formulas are used for u^ .with reference to the sign of 
d^ . The same approximations can be used for the derivatives ap- 
pearing in and 0 ^ . he recall that they are affected by coef- 
ficients which should be small; therefore, any approximation 
should serve the purpose. For the same reason, any approximation 
to P^ and appearing in the grid corrective terms in and C* 
should work. One can cnoose to use centerea differences or a 
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Extension to mul t id imen s ion* 1 problems. Unsteady flows 


MacCormacx-type of alternating backward and forward differences. 

X J 

Once P_ and P have been evaluated, P may be obtained as 
T T T 

P t «pJ ♦ Pj (36) 

Obviously, the same result is obtained by coding in one single 
equation : 


XX y v 

P ♦ A X P >A IP ♦ b X P -B X P ♦ 
T 11X1 22X2 1 1 Y1 2 2 Y2 


X _ X X 

* “u'Vxa-Vxi 1 


Y Y Y 

D (X v -X v ) ♦ C *0 
12 2 Y2 1 Y 1 1 


(37) 


This shows that the scheme does not belong to the class of frac- 
tional step schemes [12,131* The Initial separation of the equa- 
tions was introduced to allow two independent sets of charac- 
teristics to be defined, one in the (X, T) plane and the other in 
the (Y, T) plane. Attempts of the late fifties and early sixties 
to define ana use an optimal set of three-dimensional charac- 
teristics proved to be cumbersome, difficult to apply and, at 
times, unstable or inaccurate. Tne present approach may not be 
appealing from a formal viewpoint but it fulfills the purpose of 
defining the domain of dependence of the point to be evaluated as 
far as pressure and any other variable having the same domain of 
dependence as the pressure are concerned. 


I* Ifar.fig-dlafinalflnal , st eady , supersonic flows 


Three-dimensional, steady, supersonic flows must be 
treated in a different way, although the basic concepts remain 
unaltered. This is because the physical nature of the three in- 
dependent variables is the same, despite one playing a time-like 
role. We begin by observing that, after introducing two angular 
variables, a and n, defining the direction of the velocity vec- 
tor, the Euler equations can be written in the matrix form: 

f T ♦ A f x ♦ Bf y ♦ c - kP T *0 (36) 
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where 


» A « 


1 1 

*12 

a 

13 


*11 

b 12 

21 

*22 

0 

. B a 

°21 

b 22 

31 

0 

1 

CM 

CM 

m 


6 3 1 

0 


22 


C a 




.4 ■ 


We can split the system (3b) into two systems: 


k 
L 'J 


(39) 


P T * *u\ * *12°X 


o ♦ a P ♦ a a 
T 21 X 22 X 


*t 3 \ * °t 


a 0 


*3. P X 


* V Vt • 0 

♦ *22"x * V k 3 P i * 0 


C«0) 


and 



p ; 

♦ 

6 ,i p x 

♦ 

b o 
12 Y 

♦ 

b 13"x 

a 0 



Y 

°T 

♦ 

b 2, p x 

♦ 

b a 
22 Y 


- -A 

2 T 

a 0 

(41 ) 


n Y 

♦ 

6 J, P x 



♦ 

6 22V k 3 P T 

a 0 


so that, if each 

equation 

(40) is 

added 

to its 

counterpart 

in 

(41 ) , (38) 

is 

obtained again 

, provided 

that 



P ( 

s 

r 

p J 

• p J. 

°T 

X 

* 

T 

Y 

♦ 

T 

X 

* n_ * n_ 

T T 

X 

♦ 

T 

(42) 


For each new system two values of X can be found, defined by 
X 2 X 

(X ) -(a ♦a +a k ♦a k )x ♦ ( a a -a a -a a ) = 0 ( 4 3 ) 

i 11 22 133 122 i 1 1 22 1331 122". 

and 

y 2 y 

(x ) -(b *b +b k ♦b k )X +(b b -b b -b b )=0 (44) 

i 11 22 13 3 12 2 i 11 22 13 31 12 21 

respectively, and a third value, defined by 


. X 

X : a 

3 22 


(45) 
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and 


x Y .b 

3 22 


(46) 


respectively. The direction defined by X_ plays the role of a 
streaaline in (X,T) plane or the (Y,T) plane. Non-trivial sets 
of aultipllers can be found in connection with each of the X^ to 
build two sets of coapatibility equations. Specifically, 


X , X 

u s a -X 
11 22 1 


u ■ -a 

12 12 


v » - a 

13 13 


u « b -X 
11 22 1 


u « -b 

12 12 


u * - b 

13 13 


(47 ) 


for i*1,2, and 
X 


" 3 i*° 

* V22 1 

x 

u ia ♦ k a 

33 21 2 22 


/ .0 

31 

“ 3 2*' l6 31 * Vm' 

Y 

u >b ♦kb 
33 21 2 22 


(4d 1 


The compatibility equations, ir. which we write f^ and f in- 
stead of f^ and f^ (see Section 2), are 

I X . X ,X. % . X x X . X X . , X 

1 a (P ♦ X P )-a (o ♦ x o )-a (n ♦ X n ) ♦ Eu c *0 

T 1X1 12 T 1X1 13 T 1X1 1 J J 

X, X , X„ % , X X . , X , X , „ X . . 

-a (P ♦ X P )-a (o ♦ X o )-a (n ♦ X n ) ♦ lu c sO (49) 

T 2 X2 12 T 2 X2 13 T 2 X2 2J J 

,X, X „ t X , X % X , X , , X 

-4 (P ♦a^.P v ,) + u,_(o ♦a o )♦« (n ♦ a„n ) ♦ Iu„ c *0 
T 22 X3 32 T 22 X3 33 T 22 X3 3J J 


where 


a * u - u k-u k 
11 12 2 13 3 


4 ”21 k 3‘ a 31 k 2 


Similar equations can be written with respect to Y. 


(50) 
Prom the 
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first two equations (49) we get one equation for P ^ : 

X X X X X XX 

2a P ♦ a (X P ♦ X p )-a (X o -X a )- 
T 1 XI 2 X2 12 1 HI 2 X2 

* “X>v 

x 

and another equation for the unknown a : 

, * X, X , X Q » ,,X , x 

-j, . . <VirVx2 ) -*i*‘Vxi * Vxi 1 ' 

■*u“?"xi * ‘‘"xj 1 * tU M * 

where 

XX X 

• "h't * '.3 n I 

X 

Finally, considering P.^ as defined by 151), the third of 
yields 


X x X ,„ X 
* * 4 ( P_ ♦ a 


22 X3 


) ♦ »*-.• 


32 22 X 3 


X 

♦ u _ . a 


33 22 X3 


♦ «. }J 0 .0 


where 


XXX 
e » u a 
32 T 


X X 
u _ _ n_ 
33 T 


X X 

and o_, n ^ can be obtained by solving the system ( 53 ), ( 55 ). 
should proceed in a similar way for the second set (41) and 
apply (42) to get P^ a ^ and f> T . 

The rules for discretization are the same as outlined 
Section 3. The f ^ and f^ approximations coincide with f^ 

f if the signs of^X, and 5 X coincide; otherwise, they* 1 

coincide with f „ and f . 

X2 X2 


In the simple case of Cartesian coordinates with 


(51 ) 

(52) 

(53) 
(49) 

(54) 

(55) 

One 

then 

in 

and 

will 

no 
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s t retoning , 




a .*,,■0 

• * i d 



j , 1 

*2.** 

’*22’*' 

•23.°. Vj,« 

( 56 ) 

1 

'* 31 * 0 - 

* 52 '° 

* * 35 ** * ** 3**21 * 



t> n ■<!/« , 

b 12*° 

• b 13 **"* 


( 

°2 1 -0 ' 

6 22*" 

’ * 23 *°* ** 2 ** 2 l” 

( 57 ) 


V*2l’ 

,32.° 

’ 6 33 *'’’ ** 3 ** 2 l" 


V 

0. 

Then , 




X 

X s 

4 

(oS8)/k , 

x a! 

* * , 

,2. 2, 2 ,1/2 
(1 ♦ 0 )/a -1 ] 



* 


w 


( 56 ) 

' 

' 1 
X 3 
4 

( n*8 * )/« , 

x a! 

* • , 

2 2 2 1/2 
[w <1 ♦ n )/•-!] 


* 


c 

w 




i x 

*n*-* 

X 

1 

Y Y 

u s n - 1 
il i 



* n 

“13" 0 


i „ 
“l2*° 

X , 

"is**’" 


(59) 


for is 1 ,2 , and 


x „ 

“ 3 i*° 


“32“*2t” 


A - t « > 

55 1 *21° * ” 


Y _ 

M 31*° 

Y , 2 

U 32- a 21 (n 


U 33 * a 2l n ° 


1 ) 


(60) 


so that 
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7e 


X, 

t« t /« 


Y 

♦ « tn 


Y 

Y 
T 


/« 

/« 


X a f ,2 , . X X. Y i 

* * ,l(o ♦ 1 ) n^-o no^ J , * ■ 


yw 


YW 


-r * , 2 
^lonn^-Cn 


Do*] 


( 61 ) 

( 02 ) 


(03) 


All comp 1 ications appearing in the general formulation 
above depend on the geometry of the computational grid and should 
not be relevant from a physical viewpoint. 


A. laaifl aafl reauita 


Some examples should make clear the ability of the scheme 


to represent difficult transonic flows. he will begin with two 
one-dimensional cases, one of which is -'nsteady and the other 
steady and supersonic. 



Fig. 7 

The first example is illustrated in Fig. b, wnich snows 
the Mach number distribution over a stretch of an infinite duct 
having a variable cross-section (Fig. 7). In the portion of the 
pipe considered in the figure the flow, which started from rest 
and passed through an unsteady evolution, has already reached a 
steady state. A shock was originated and moved to a final steady 
position. The calculation used to get Fig. 8 has no provision 
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for fitting a shock; therefore the flow is homen tropic . 
Nevertheless, the results are surprisingly close to the results 
obtained using a shock-fitting technique (Fig. 9) which, in turn, 
coincide with the close-form solution of a steady, quasi-one- 
dimensional flow. A single figure cannot provide too many de- 
tails; a more minute analysis can be found in [10] where the en- 
tire range of Mach numbers is explored and the unsteady evolution 

leading to the asymptotic steady state ir analyzed. I would like 
to warn the reader against concluding, from a glance to Fig. 6, 

that the l-scheme has shock-capturing virtues; such a feat is 
made impossible by the inability of the scheme to provide entropy 
Jumps. It is undoubtedly true, however, that the scheme remains 
stable and accurate under very exacting transonic circumstances, 
even when the computational mesh is very coarse; in Fig. 10 
asymptotic results obtained with the mesh of Fig. 6 and a much 
coarser mesh are shown, and no deterioration can be observed. 

The second example is illustrated in Fig. 11, which shows 
the Mach number distribution in a channel opening into vacuum. 



Fig. 9 


Tne flow at the channel inlet is unifora; its Mach number equals 
1.04. The calculation [14] is performed advancing on a grid of 
which Fig. 12 shows typical l*nes. On each grid line crossing 
the channel, 21 points are considered. In this case, the accura- 
cy of the scheme and its ability to handle difficult situations 
is deaonstrated by the perfect splitting of the coaputed flow 
into two regions, a siaple wave adjacent to the lower wall and a 

source-like pattern in the re of the channel. Such a splitting 
is due to the fact that no down-running characteristic can reach 

the lower wall. In Fig. 11, line AL is the first significant 
down-running characteristic. 


An impressive example of the ability of the scheme to 
handle two-dimensional , unsteady flows is shown in Fig. 13. In 
it, we see the steady shock-and-isobars pattern reached asymptot- 
ically about an ablated blunt body via an unsteady calculation. 
The free stream Mach number is 12. The computational grid 
(wrapped around the body by conformal mapping [15]) has 30 inter- 
vals along the body and 15 intervals between body and shock. In 
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Fig. 10 


the figure, in addition to isobars, two sonic lines are shown. A 
thin subsonic layer appears in the concavity of the body. The 
bow snock shape matcnes the one obtained experimentally [loj. 

Finally, an example is given of the application of the 
technique to three-dimensional, steady, supersonic flows. In 
Fig. 14 the bow shock and the isobar pattern about an elliptic 

cone, with a 14:1 axis ratio and a leading edge angle of 20 are 

shown. The free stream Mach number is 2 and the angle of attack 
o 

is equal to 5 . Despite the coarseness of the computational mesh 
(24 intervals around the body ana 12 between body and bow shock) 
the pressure distribution on the body is remarkably accurate, as 
shown in Fig. 15. The results are part of an extensive series of 
tests involving vehicle shapes of increasing complexity, to be 
published in the near future. 
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